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THE EIGENVALUE SPECTRUM OF THE ORR-SOMMERFELD PROBLEM 


I. INTRODUCTION 


This report is primarily concerned with a numerical investigation of the temporal 
eigenvalue spectrum of the Orr-Sommerfeld equation. The solution and the general 
properties of this equation are still a topic of interest nearly seven decades after its 
formulation. The long-standing interest in the problem is primarily due to two reasons. 
The first reason is due to physical considerations in that the equation has served as a 
successful model for the study of the linear stability of one class of hydrodynamic flows, 
namely, that of parallel flows. The second, and more abstract, reason is due to the 
complex mathematical structure of the equation, whereby this equation, together with 
the one-dimensional Schroedinger equation, belongs to a class of linear differential 
equations that contain a large parameter and turning points. This equation is further 
complicated by the fact that it is a non-self-adjoint problem. The derivation of the 
Orr-Sommerfeld equation for the linear stability problem together with an excellent 
discussion of its solution can be found in the classical monograph by Lin [ 1 ] . A 
comprehensive historical survey on the solution and the properties of a whole class of 
differential equations to which this problem belongs is discussed by McHugh [2]. 

In this report only one property of this problem which has received attention 
recently will be studied, namely, the type of the eigenvalue spectrum of the equation. The 
knowledge of whether this equation possesses an infinite or a finite spectrum is necessary 
for one method of study of the nonlinear hydrodynamic stability problem because, as 
shown by Eckhaus [3], the solution of the nonlinear problem may be obtained through 
an expansion in terms of the eigenfunctions of the linear problem. Hence, it is necessary 
to show that the linear problem possesses a complete set of eigenfunctions. Also, there is 
the more fundamental issue which requires the proof that the eigenfunctions of the linear 
problem form a complete set. This is necessary to determine the class of initial 
disturbances that may be investigated by the method of normal modes. 

In principal, the completeness property has been neglected in the past for many 
reasons, .and one of them may have been that a completeness proof for the general 
non-self-adjoint problems simply does not exist in the mathematical literature. However, 
the first to brave that hurdle was Schensted [4] , who proved an expansion theorem for 
the Orr-Sommerfeld equation for two specific flow profiles, namely, the plane Poiseuille 
flow and the pipe flow. This effort was followed sometime later by DiPrima and Habetler 
[5] who proved a completeness theorem for the general linear stability problem for 
bounded flow profiles. However, to date there is no proof of a similar theorem for the 
stability problem of the boundary layer flow because this flow profile is semi-infinite; 
i.e., the flow is bounded by a wall on one side and extends to infinity on the other. 
Also, the controversy surrounding the nature of the eigenvalue spectrum of the 
Orr-Sommerfeld equation for the boundary layer flow has not been settled yet. 



Of course, the only way to settle this issue is by a rigorous mathematical proof 
similar to the one given by Schensted for the bounded flow situation. However, it is 
possible to explore this question, in principle at least, through an investigation similar to 
the one presented herein. Indeed, two such attempts have been undertaken recently, one 
by Jordinson [6] and the other by Mack [7]. In these studies it was shown through a 
search in a finite portion of the eigenvalue plane that there exist only a few eigenvalues 
of the equation, leading to the conjecture that the eigenvalue spectrum for the boundary 
layer profile is finite. However, the number and location of these eigenvalues varied from 
one search to the other. In this report it will be shown that the eigenvalue spectrum of 
the Orr-Sommerfeld equation for the boundary layer profile may be infinite. Also, this 
investigation will show that the findings of both Jordinson and Mack are not 
contradictory within the framework of these respective studies. 

Because of the complexity of the equation, closed form analytical solutions are 
precluded in this study and the problem will be investigated through numerical 
approximations only. There are many numerical techniques in the literature for solving 
the Orr-Sommerfeld eigenvalue problem, but none is suited for the present investigation. 
Most of these methods fall within two general classes, the matrix method and the 
shooting method. In the former, the differential eigenvalue problem is reduced to an 
algebraic problem through various approximations to the solution. The most widely used 
of these methods is the finite difference approximation technique such as the one used 
by Jordinson [8]. Another example of this class is the Chebyshev Polynomial 
Approximation method such as the one employed by Benek [9]. However, these 
methods possess two shortcomings that preclude their use in the present investigation. 
One is that the number of eigenvalues sought and determined is directly proportional to 
the number of mesh points used or the number of terms retained in the polynomial 
approximation. The second deficiency is that these methods do not provide for an a 
priori control on the section of the complex eigenvalue plane that needs to be searched. 

As for the shooting method, the boundary value problem is solved as a system of 
initial value problems by choosing a suitable guess on the missing initial conditions. The 
equations are then numerically integrated to the desired end point, where the guessed 
initial conditions are corrected through the requirement that the equations satisfy the 
given end condition. For this purpose, an iterative procedure such as the 
Newton-Raphson-Kantorovich method may be used. In linear problems such as this one, 
only one iteration is necessary to determine the missing initial conditions and eigenvalues. 
This method has been extensively used by many investigators, such as Mack [10] and 
Davey [11]. However, the technique possesses a major shortcoming in the sense that only 
one eigenvalue may be determined at a time, and in order to start the iterative process 
and have the iterations converge, it is necessary to have a good estimate of the 
eigenvalues that need to be determined. This is obvious since most iterative processes are 
local in nature. Previously, such good estimates on the eigenvalues were determined from 


2 



the known approximate analytical solutions or by a variational method such as that of 
Lee and Reynolds [12]. However, when' such estimates do not exist or are cumbersome 
to produce, the method will involve trial and error procedures that may be very costly 
and frustrating. 

It is clear from the literature survey that the known numerical techniques are not 
adequately suited for the present investigation and that a new method must be developed 
to carry out the desired study satisfactorily. Such a method is developed here. It belongs 
to the shooting method class and has the additional advantage that it is capable of 
locating more than one eigenvalue at a time and does not require an initial guess. The 
method is global in the sense that all the eigenvalues that lie within a finite region of the 
eigenvalue plane are identified and determined. The method is also automatic in the sense 
that once the boundaries of the domain desired to be searched are given, it proceeds to 
locate the eigenvalues without any further information from the operator. Also, the 
technique has the built-in capability of determining multiplicities of the eigenvalues if any 
exist. 


II. THEORETICAL BACKGROUND 


In this section the theoretical background that is necessary for the development 
of the numerical technique for the present investigation is outlined for the general linear 
differential eigenvalue problem. However, before entering into this discussion, the specific 
problem that will be studied later will be defined so that it may serve as a guide in the 
development of the technique. 

The specific equation that will be studied is the one governing the amplitude 
function of a two-dimensional, wave-like, infinitesimal perturbation in an otherwise 
stationary, two-dimensional flow, namely, 


(D 2 - a 2 ) 2 (ft = i cf R [(u - c) (D 2 - a 2 ) (ft - (ft D 2 u] 


( 1 ) 


where D denotes differentiation with respect to the spatial coordinate y. Equation (1) is 
commonly known as the Orr-Sommerfeld equation in which (ft is the complex amplitude 
function; a and c are the disturbance wave number and speed, respectively, while u(y) is 
the stationary velocity profile that is assumed to be defined throughout the range of 
integration. R is the Reynolds number which is defined by UL/p, where v is the 
kinematic viscosity. All of the variables in equation (1) have been made dimensionless 
with respect to a suitable length scale L and a velocity scale U. 

The appropriate boundary conditions for equation (1) which are physically viable 
are either the no-slip conditions at a solid boundary or the requirement that the 
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perturbation decay to zero at infinity when the fluid is infinite in extent. Such boundary 
conditions are commonly known as homogeneous,' and the specific conditions for the 
cases that will be studied will be specified separately later. Equation (1) and the 
homogeneous boundary conditions form an eigenvalue problem, and it is our task herein 
to investigate the nature of the eigenvalue spectrum which belongs to this problem. 

Equation (1) may be written in a general operator form, 


L! 0 + c L 2 0 = 0 , (2) 

where 

Lj = (D 2 - a 2 ) 2 - i a R [u(D 2 - a 2 ) - D 2 u] (3a) 

L 2 = iaR(D 2 - a 2 ) , (3b) 

where the wave speed c = c r + icj is taken to be the eigenvalue. However, any other 
combination of the parameters that appear in equation (1) may be taken as eigenvalues. 

Let us consider a general differential eigenvalue problem of the form (2); i.e., 

Li 0 + X L 2 0 = 0 , (4) 


where L! and L 2 are linear differential operators of the form, 


L(0) = p 0 (y) D n 0 + p 1 (y)D n “ 1 0 + •■■ + p n (y) 0 


(5) 


For the cases of interest here, we need to limit outselves for the case where the order of 
L 2 is less than or equal the order of Lj . X in equation (4) is the complex eigenvalue, and 
0 is the complex eigenfunction of the real variable y. Let us assume that the coefficients 
of equation (5), pj^y), possess continuous derivatives up to order n on the interval of 

integration [a,b]. Also let the differential equation (4) be subject to the following 
homogeneous boundary conditions: 


Bj(0) = 0 (i n) 


( 6 ) 
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which are specified at both ends of the interval [a,b] . Note that the boundary conditions 
Bj(0) may be composed of linear combinations of the 2n quantities. 


0(a), D0(a), • • 

• . D n 'h(a) 

0(b), D0(b), • • 

• . D n_1 0(b) 


Let us now define the determinant f(X) by 


f(A) = det 


Bi (00 
B n (0 ! ) 


Bi (0 n ) 
B n »n> 


(7) 


where 0j,02, . . . , 0 n is a- fundamental system of Cauchy solutions of equation (4) 

collectively satisfying linearly independent sets of initial conditions at y = a. If the 
coefficients of the differential operators L x and L 2 , namely, p^fy), and the boundary 

conditions Bj(0) are independent of the eigenvalue A, the following theorem concerning 

the eigenvalues of equation (4) may be invoked, Naimark [13]: The eigenvalues of the 
operators of equation (4) together with the boundary conditions of (6) are the zeroes of 
the function ff\). If /f A) vanishes identically, then any number A is an eigenvalue of the 
problem defined by (4) and (6). If, however, f(\) is not identically zero, equation (4) has 
at most denumerably many eigenvalues and the eigenvalues have no finite limit point. 

It can also be shown [13] that when the differential operators L 1 and L 2 are not 
singular, then f(A) is an integral analytic function of A. This last property is of primary 
importance for the development of our numerical scheme, as will be shown. 

Thus, at least in principle, the problem of locating the eigenvalues of equation (4) 
together with the boundary conditions given by (6) is determined. However, determining 
the zeroes of the function f(A) is more complicated, and the common procedure in the 
past has been to determine only one zero at a time. This was accomplished, as explained 
in Section I, through local iterative techniques that require a good initial guess to begin 
them. However, if one is able to generate a global technique for determining the zeroes 
of f(A), one possesses the capability of determining more than one, or all, of the 
eigenvalues of the problem defined by (4) and (6) at one time. Toward this end one may 
use the following theorem from the theory of complex analysis which provides the means 
for locating all of the zeroes of a complex function in a finite region of the complex 
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plane. This theorem may be stated in the following way: If y is a regular curve in a 
simply connected open set not passing through any zero or pole of a meromorphic 
function f(\) in that set, the integral 


1 

2m 


f 


f(X) 

W 


dX 


( 9 ) 


is equal to the excess of the number of zeros over the number of poles of f(\) inside y (a 
pole or zero of order k counting as k poles or zeros, respectively ). 

Since, as mentioned earlier, f(X) in our specific problem is an analytic function of 
X, then within any closed region in the X-plane there does not exist any pole; hence, 
expression (9) may be modified to read 


__L 

2m 


\ m 


dX 


N 


( 10 ) 


where N is the number of zeros within the curve y. Furthermore, the above theorem 
admits a simple generalization that may be valuable later. Again, let f(X) be analytic 
inside and on y and let co(X) be also analytic inside and on y. Then we have the 
following result: 


s / Sr dx = j; qr «<v • 

where the sum is taken over the zeros of f(X), and q r is the order of the zeros X f . The 
proof of expression (11) may be found in Goodstein [14]. 

To carry out our investigation, the previous two theorems together with 
expression (11) will be used in the form of numerical approximations. The details of such 
approximations will be discussed in Section III. 
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IM. NUMERICAL DETAILS 


A. The Numerical Evaluation of the Zeros of a Complex Function 

As outlined in Section II, the procedure for locating the eigenvalues of the 
differential eigenvalue problem reduces to locating the zeros of the complex function f(c) 
provided that f(c) could be determined at an arbitrary number of points {hereinafter c 
will replace X used previously). The procedure for the numerical evaluation of f(c) and 
f'(c) for any value of c will be discussed later. Note that because of the complexity of 
equation ( 1 ), f(c) can only be evaluated numerically; however, for a less complex 
equation it might be possible to evaluate this function in a closed form and thus alleviate 
the necessity for a numerical investigation. Assuming temporarily that f(c) andf'(c)may 
be evaluated around any contour 7 , then the zeros of that function may be located 
through the use of expressions ( 10 ) and ( 11 ) in the following manner (a detailed 
discussion on the technique that is outlined here may be found in Reference 15). 

First, the number of zeros inside a prescribed region of the complex c-plane, here 
this region is taken to be a circle, may be determined numerically through an appropriate 
quadrature formula approximation for expression (10). Once the number of zeros inside 
7 , N for example, is determined, their location may be found in the second step. By 
letting the function te(X) in ( 11 ) take successively the values c, c 2 ,c 3 , . . . , c^, all the 
power sums of the zeros inside 7 may be determined by performing the following N 
contour integrations: 


r f'(c) . v 

1 = J c fW Zj c j 

7 IW j=l 


N 


s, = / c * IM dc = 2 

7 j=l 

fU \ N 

s n = f cN Sr dc = 2 

7 j=l 


N 


( 12 ) 


N 


where cj denotes all the zeros of f(c) inside y. Again, the numerical values of the contour 
integrals in ( 12 ) may be obtained through an appropriate quadrature formula. 
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Using the complex numbers Sj S2> • • . , s^, a polynomial P(c) of degree N may be 

constructed whose zeros inside y coincide with the zeros of f(c). Note that this 
polynomial is not an approximation of the function f(c). For the construction of such a 
polynomial, Newton’s formula is used in which the coefficients of the polynomial P(c) 
are evaluated from the power sums of the roots of the polynomial. Once the equivalent 
polynomial, P(c), is constructed, any polynomial root finder may be employed to extract 
the zeros of P(c). The zeros of f(c) are coincident with the roots of P(c) inside y. The 
specific polynomial root finder subroutine used in the present investigation utilizes the 
quadratic method. It should be mentioned that the numerical technique for finding the 
zeros of the function f(c) as described previously is also capable of determining 
multiplicities of the zeros if any exist. However, multiple zeros were not found in the 
present problem within the domain of search. 


As for the numerical evaluations of the quadrature integrals in (10) and (12), the 
following formula developed by Delves and Lyness [15] for the integration around a 
circle in the complex plane was employed. Let the contour of integration, y, be the circle 
7(c 0 ,r) whose center is c 0 and radius is r; upon translating the center of the circle to the 
origin, these quadrature expressions may be written as integrals with respect to the 
parameter t in the following way: 


S N 


i 

f exp [27ri(N + l)t] r 
0 


f' [cq + r exp (27rit)3 
f[cQ + r exp (27rit)] 


(13) 


where N, which is the number of zeros inside y , is set equal to zero when evaluating the 
quadrature integral of (10). Upon using the trapezoidal rule for approximating the 
quadrature integral in (13), we obtain the following numerical formula for the integration 
around the circle 


S N 


1 m 

— J] '/'n (i/m) + 0 ( Am ) 
m j=l 


(14) 


is the integrand in (13), and m is the number of discrete points around the circle y 
at which f'(c)/f(c) is evaluated. A in (14) is defined by 


A = max (A 1 , A 2 ) 


(15) 
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where 


|A| < 1 
A j = max v^ x 


~ max r/r j 


where rj is the distance between the center and the closest zero to the perimeter of the 

circle 7. Formula (14) is convenient in that by increasing the number of function 
evaluations around the perimeter of y from m to 2m requires only m additional function 
evaluations. Also, it is notable that the convergence of (14) to (13) is linear; i.e., for each 
additional function evaluation the error is reduced by a constant factor. 


B. The Numerical Evaluation of f'{c) 

An examination of the quadrature expressions (10) and (11) reveals that the 
evaluation of the derivative of f(c) with respect to c is necessary for the implementation 
of the root finding procedure. This derivative may be obtained by the usual rule of 
determinant differentiation that requires the evaluation of the derivative of each element 
of (8) with respect to c. To achieve this, the differential system (2) is augmented with 
yet another system of the same order governed by the following differential equation: 


Li + c L 2 4* + L 2 0 — 0 


(16) 


where <h = 30/dc. The differential operators and L 2 in equation (16) are the same as 
those defined by (3). Note that since f is an analytic function of c, the differentiation 
may be taken with respect to either c r or q without any loss of generality. Equation (16) 

is then solved subject to the following initial conditions: 


<h(a) = D<h(a) = 0 = D 2 <f>(a) = D 3 4>(a) 


(17) 
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C. The Initial Value Integrators 


It is obvious from the previous discussion that to determine the eigenvalues to 
equation (1) inside the contour y, the functions f(c) and f'(c) need to be determined 
numerically at discrete points along the perimeter of the circle y. These values may be 
obtained by evaluating the determinant in expression (8) for f(c) and a similar 
determinant for f'(c), where each column of the determinant represents a Cauchy 
solution of equation (1). Naturally, each Cauchy solution involves an initial value 
integration of equation (1) over the interval [a,b] with the appropriate initial condition. 

The forward integrators used in the present study are of the 

Runge-Kutta-Nystrom type for a system of second-order ordinary differential equations. 
This integrator utilizes the fact that the differential operators and L 2 do not possess 
any odd order derivatives, leading to a reduction of the computation time of the whole 
system by approximately one-half over the classical Runge-Kutta methods for first-order 
equations. While many extensive higher order formulas with step size control exist [16] , 
the specific integrator used here is a seventh-order Runge-Kutta-Nystrom with a fixed 
step size. 

The Orr-Sommerfeld equation (1) presents difficulty when it is numerically 
integrated as an initial value problem because some of the parameters that appear in the 
equation are widely separated. This is apparent when it is observed that the wave number 
a is of order 10" 1 while the Reynolds number for any problem of interest is of order 10 3 
or more. Since in the numerical integration of this equation one begins with a system of 
linearly independent (orthogonal) initial vectors, because of the wide separation of the 
parameters, these vectors become increasingly parallel as the integration proceeds to the 
final point, resulting in their being linearly dependent and rendering the matrix for the 
determinant (8) ill-conditioned. 

Few methods have been devised to overcome this situation, and, to the author’s 
knowledge, they are all implementations of one basic idea. Since the initial vectors were 
orthogonal, it is necessary to ensure that the solution vectors remain orthogonal 
throughout the range of integration. This is achieved by reorthonormalizing the solution 
vectors at discrete points along the path • of integration. Such a technique was first 
developed by Godunov [17], which was later improved by Conte [18], whereby the 
Gramm-Schmidt orthonormalization technique is employed. In the present study, 
however, because equation (16) is also being integrated simultaneously with (1), one must 
ensure that the solution vectors of (16) are transformed in the same manner as those of 
(1). That is, if the solution vectors of (1), for example 0 t 0 2 , • • • > 0 n , are transformed to 
> 02 5 • • • » 0 n trough the transformation 


01 


0i 

02 

= [P] 

02 

$ 

_ 


0 

_n_ 


( 18 ) 
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at some point along the integration interval the solution vector of (16), for example 
4> n , should be transformed at that point by 


*1 



*2 

= [P] 

i 2 

L n J 


$ 

L nJ 


(19) 


where [P] is the orthonormalization transformation matrix at the point in question. 


IV. THE PLANE POISEUILLE FLOW 


A. Definition of the Problem 

This type of profile is distinguished by the fact that the flow is bounded between 
two parallel stationary planes located at y = ±l. The stationary laminar profile for this 
configuration is a parabolic function of y given by 


u(y) = 1 - y 2 . (20) 

The length scale L for this problem is half the channel width, while the velocity scale U 
is the maximum velocity midway between the planes. The boundary conditions 
appropriate to equation (1) for this case are the usual no-slip conditions at the walls, 
which may be written as 


0(± 1) = 0 = D0(±1) 


( 21 ) 


Since the mean flow profile in expression (20) and the boundary condition are 
symmetric in y, then the solution for equation (1) may be split into odd and even parts. 
Thus, for solving the present eigenvalue problem it is sufficient to consider either the odd 
or even solutions. Splitting the problem in this manner permits the solution of equation 
(1) for the half-range [0, 1] in place of the full range [-1, 1] by replacing the condition 
at one wall with the following conditions at the center of the channel: 
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D0(O) = D 3 0(O) 


for the evert solution 


(22a) 


and 


0(0) = D 2 0(O) 


for the odd solution 


(22b) 


thus permitting the reduction of the integration interval by one-half. Although the 
present technique of eigenvalue search applies equally well for either the even or the odd 
solution, only the even solution will be investigated here. 

With the problem and boundary conditions thus defined, the determinant in 
expression (8) will then take the following form: 


DM0) 

002 (0) 

D0 3 (O) 

D0 4 (O)' 

O 3 0, (0) 

D 3 0 2 (O) 

D 3 0 3 (O) 

D 3 0 4 (O) 

0,0) 

0 2 O) 

03 (1) 

04(1) 

00,(1) 

D0 a (l) 

D0 3 ( 1 ) 

O0 4 (1) J 


(23) 


This determinant is of high order requiring the solution of four Cauchy problems and 
may be reduced considerably in the following manner. Since the boundary conditions 
(21) and (22) are separable with two of them applying at y = 1 and the remaining two at 
y = 0, a reduction of the order of the determinant may be achieved by solving an initial 
value problem defined by augmenting the given conditions at y = 0 by two further sets of 
two linearly independent conditions at that point. When the resulting fundamental 
solutions are 0, and 0 2 , where 0j = (0j, D0p D 2 0-, D 3 0j), expression (23) may then be 

replaced by the characteristic determinant equation expressing the compatibility of a 
nontrivial combination of these solutions and their derivatives with the given boundary 
conditions at y = 1 ; namely, 


f(c) - det 


Ml) 

00 ,( 1 ) 


MD 

DMD 


(24) 


Thus, the solution of two rather than four Cauchy problems is sufficient, resulting in a 
great labor savings. The linearly independent initial vectors for these two Cauchy 
problems are: 
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•MO) - (1,0, 0,0) 


* 2 < 0 ) = ( 0 , 0 , 1 , 0 ) 


( 25 ) 


Thus, the forward integration is started at y = 0, with (25) as initial conditions; f(c) is 
then determined by evaluating the determinant given by (24) at y = 1 . A similar 
procedure is followed for the determination of f'(c). 


B. The Search Technique 

The following search procedure in the eigenvalue plane was employed in this 
investigation. The region desired to be searched was spanned by circles of uniform radius 
as shown in Figure 1 in such a way as to cover the entire desired search area by a 
minimum number of circles. Although it is desirable to span the complete domain of 
interest by one circle, it was found that there exists an upper limit on the size of the 
individual search contour. It was observed that whenever the radius of the circle in 
question was greater than a certain value, the quadrature formula as given by (14) failed 
to converge. The value of the maximum radius allowable was found to be dependent on 
the Reynolds number in equation (1); i.e., the smaller the Reynolds number was, the 
larger was the allowable value of the radius. 

Thus, using the technique outlined previously, each circle was then searched for 
eigenvalues within it. To begin the individual search, f(c) and f'(c) were evaluated initially 
at eight equally spaced points around the circumference of the circle. The number of 
points was then doubled, and the quadrature value obtained from the 8- and the 16-point 
approximations of (14) were compared. When the difference in the value of s^ between 

the two iterations was less than 1 percent, the quadrature value was considered 
converged. If, however, this criterion was not met, the number of function evaluations 
around the perimeter was doubled and the convergence criteria were tested again. 

All of the computations reported hereafter were performed on the IBM 360/65 
computer in double precision arithmetic when possible. Since the time required for each 
function evaluation is considerable (for each such evaluation, the integration of 16 real 
second-order differential equations over the interval [0,1] was necessary), the 
computation time for the search of each circle is directly proportional to the final 
number of function evaluations around the perimeter. For example, the time required for 
128 such evaluations was approximately 3 minutes of computation time. More details on 
computation times and accuracies may be found elsewhere [19] . 
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Figure 1. A portion of the complex c-plane indicating the method of search for 
eigenvalues within a specified region of that plane (X indicates eigenvalues 
found inside the specific circles). 


C. Results 

As discussed in Section I, the eigenvalue spectrum of the Orr-Sommerfeld 
equation for the plane Poiseuille flow has been investigated in earlier works, but none of 
these investigations dealt in depth with an eigenvalue search. The only complete search to 
date is that of Mack. Employing a different search scheme from the present one, Mack 
has tabulated the 32 known symmetric eigenmodes for a~ 1.0 and R= 10 000 that lie 
within the rectangle 0<c r <1.0 and -1.0<Cj<0 in the complex c-plane. Hence, the 

results of the present scheme for that particular case will be compared with Mack’s 
tabulation. 

All of the eigenvalues that were found in the region 0 < c r < 1.0 and -1 .0 < Cj < 0 

for the parameters a= 1.0 and R= 10 000 are tabulated in Table 1 and plotted in Figure 
2. These eigenvalues agree as to location and number with those of Mack, and no 
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TABLE 1. EIGENVALUES FOR THE PLANE POISEUILLE FLOW AT a = 1.0 AND 
R = 10 000 IN THE REGION OF THE c-PLANE BOUNDED BY 
0.0 < c r < 1.0 AND -1.0 < cj < 0.0. 


A Family 

P Family 

S Family 

c r 

c i 

c r 

c i 

c r 

c i 

0.2381 

+0.0038 

0.9636 

-0.0362 

0.6776 

-0.3437 

0.3487 

-0.1245 

0.9383 

-0.0614 

0.6745 

-0.3898 

0.1894 

-0.1828 

0.9071 

-0.0922 

0.6673 

-0.4316 

0.4749 

-0.2087 

0.8798 

-0.1194 

0.6723 

-0.4833 

0.3685 

-0.2388 

0.8515 

-0.1474 

0.6716 

-0.5324 

0.5871 

-0.2672 

0.8231 

-0.1755 

0.6710 

-0.5833 

0.5129 

-0.2866 

0.7948 

-0.2035 

0.6704 

-0.6359 

0.6829 

-0.3076 

0.7665 

-0.2316 

0.6700 

-0.6903 

0.6361 

-0.3252 

0.7381 

-0.2597 

0.6720 

-0.7436 



0.7089 

-0.2877 

0.6692 

-0.8044 





0.6689 

-0.8642 





0.6687 

-0.9258 





0.6685 

-0.9893 


different eigenvalues were found in this region. Following the nomenclature adopted by 
Mack, these modes may be classified into three distinct families; the A, P, and S families. 
This classification is somewhat natural because of the peculiar shape of the eigenvalue 
spectrum. Two parts of the spectrum extend along approximately straight lines; these are 
labeled the P and S families. The remaining part of the spectrum does not conform to 
such a pattern and will comprise the A family. 

It is the members of the S family that are of interest in the speculation of 
whether the spectrum is finite or infinite. The eigenmodes of this family lie along a 
vertical line at the approximate location of c f = 0.667 and extend to as large a value of 

-Cj as one cares to pursue. Thus, the members of this family are inferred to be infinite in 

extent, which is in agreement with the theoretical conclusion that the spectrum of this 
problem is infinite and discrete. 

From the practical point of view, it is the members of the A family that are of 
interest because it was found that at least one eigenvalue of this family is always the least 
stable. This is apparent when one compares the spectra for various Reynolds numbers and 
a fixed a. 
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Figure 2. All of the eigenvalues found within the indicated boundaries in the c-plane for the 
plane Poiseuille flow at a = 1.0 and R = 10 000 (• indicates members of the A 
family; □ indicates members of the P family; and A indicates 
members of the S family). 

To form a general picture of the behavior of the eigenvalue spectra at various 
values of R and a fixed a, the eigenvalue spectra for ol= 1 and R= 10 2 , 10 3 , and 
6 X 10 3 are presented in Figures 3 through 5 for comparison. Also, these eigenvalues are 
tabulated in Table 2. It is obvious from comparing these spectra that the density of 
eigenvalues in all the families increases with the increase in Reynolds numbers. However, 
the general character of the spectrum is invariant with Reynolds number. Also, in all of 
these spectra the infinite character of the S family is preserved. 
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Cr 


Figure 3. All of the eigenvalues found within the indicated boundaries in the c-plane 
for the plane Poiseuille flow at a = 1.0 and R = 100. 

V. THE BLASIUS BOUNDARY LAYER FLOW 


A. Definition of the Problem 

Although the differential equation investigated for the Blasius boundary layer 
flow is the same as for the plane Poiseuille flow (differing only in the character of the 
coefficients of the equation), these two problems are completely distinct from the 
mathematical point of view. The differences are in the boundary conditions that need to 
be satisfied. In the plane Poiseuille flow the boundaries are at finite distances, while in 
the Blasius boundary layer flow one boundary extends to infinity. This last feature 
distinguishes most external flows from internal flows which are characterized by finite 
boundaries. 
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0.0 


- 0.2 


- 0.4 

C i 

- 0.6 


- 0.8 


- 1.0 


Figure 4. All of the eigenvalues found within the indicated boundaries in the c-plane 
for the plane Poiseuille flow at o: = 1.0 and R = 1000. 

Another difference in the two flows that needs to be emphasized concerns the 
parallel flow assumption that was used to derive equation (1). While this assumption is 
appropriate for the plane Poiseuille flow and a host of other physically realizable flows, it 
is not a strictly valid assumption for the boundary layer flow. However, results derived 
from equation (1) under this assumption for the boundary layer flow have been used in 
the past as a close approximation to the real problem. It is with this reservation that we 
approach the analysis for the Blasius boundary layer profile. 

In the past many results have been built on the premise that equation (1) 
possesses an infinite set of eigenfunctions with their respective eigenvalues for the Blasius 
boundary layer flow. Hence, it is considered both instructive and worthwhile to show 
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c r 


Figure 5. All of the eigenvalues found within the indicated boundaries in the c-plane 
for the plane Poiseuille flow at a = 1.0 and R = 6000. 

whether such a spectrum exists. As far as the present analysis is concerned, the 
distinguishing feature of the problem lies in the fact that the flow is only bounded at one 
wall and is free at the other end. A flow of this type is commonly known as 
semi-infinite. The mean velocity profile for the Blasius boundary layer flow is given by 

u(y) = F'(y) , (26) 

where F(y) is the stream function that is given by the well-known Blasius equation 

F'"(y) + K 2 F(y) F"(y) = 0 (27) 
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TABLE 2. EIGENVALUES FOR THE PLANE POISEUILLE FLOW AT a = 1.0 AND 
VARIOUS R VALUES IN THE REGION OF THE c-PLANE BOUNDED BY 
0.0 < c r < 1.0 AND -5.0 < Ci < 0.0 FOR R = 100, AND 0.0 < c f < 1.0 

AND -1.0 < c { < 0 FOR R = 1000 AND 6000. 




a 

= 1.0 



R 

= 100 

R = 

= 1000 

R = 

= 6000 

c r 

c i 

c r 

c i 

c r 

c i 

0.47849 

-0.16295 

0.34628 

-0.04213 

0.25980 

+0.00032 

0.69316 

-0.38258 

0.42700 

-0.29538 

0.91782 

-0.08151 

0.66166 

-0.87322 

0.64309 

-0.27822 

0.88128 

-0.11765 

0.66351 

-1.56515 

0.70554 

-0.29722 

0.40027 

-0.15257 

0.66491 

-2.45548 

0.67397 

-0.44457 

0.84473 

-0.15378 

0.66556 

-3.54274 

0.79849 

-0.19757 

0.80817 

-0.18991 

0.66589 

-4.82693 

0.88813 

-0.11008 

0.22543 

-0.20621 



0.66924 

-0.59967 

0.77159 

-0.22605 



0.66805 

-0.77344 

0.54890 

-0.24377 



0.66755 

-0.96552 

0.73482 

-0.26242 





0.43985 

-0.26845 





0.89266 

-0.29762 





0.67834 

-0.31461 





0.60960 

-0.32006 





0.67598 

-0.37438 





0.67322 

-0.43342 





0.67199 

-0.49486 





0.67110 

-0.55915 





0.67039 

-0.62632 





0.66980 

-0.69643 





0.66933 

-0.76953 





0.66893 

-0.84566 






0.66860 

-0.92484 
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Equation (27) has been made dimensionless here in such a manner as to let it 
have the same independent variable, y, as the one in the perturbation equation (1). The 
parameter K in equation (27) relates the well-known Blasius similarity variable rj to y 

by r? = Ky, where K = J^lol r? - F(r?)3 . The reference length for this case is the boundary layer 

displacement thickness 6 lf and the reference velocity is the free stream velocity. Hence 
the Reynolds number is the one based on the boundary layer displacement thickness. 

The boundary conditions appropriate for this problem are the usual no-slip, 
no-penetration conditions at the wall given by 


0(0) = 0 = D0(O) 


(2j3) 


together with the requirement that the perturbation function be bounded for large values 
of y; i.e., 


0(°°) , D0(«>) -> 0 as y -»■ °° 


(29) 


Upon inspecting conditions (28) and (29), it appears that it is desirable to 
evaluate the secular determinant (8) in terms of the wall conditions only. Here an 
argument for the reduction of the order of the determinant from four to two may be 
followed in a manner analogous to the one presented for the plane Poiseuille flow 
problem of Section IV. To evaluate the secular determinant at the wall, the forward 
integration of the two Cauchy problems is begun at the free stream end, where the 
asymptotic properties of equation (1) are exploited and carried on to the wall end where 
conditions (28) are valid. The free stream initial conditions may be generated by 
observing that equation (1) is extremely simplified for large y; namely, as y -* ®°, u(y) -> 1 
and D 0 2 u(y) -> 0, leading to a constant coefficient equation whose solution may be 
written explicitly in terms of the following four exponential functions: 


0i 2 (Yi) - exp (± oyO 

3 

03 , 4 (yi) = exp(±pyj) , (30) 


where 


p = [a 2 + i a R (1 - c)] 1/2 , Real (p) > 0 
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In (30) y x is some arbitrarily chosen large value of y that will be discussed 
further in the next subsection. Now, condition (29) may be used to eliminate the two 
solutions with the positive exponential sign in (30). Thus, the two Cauchy problems may 
be integrated from the free stream starting with the following two linearly independent 
initial vectors: 

01 (y i ) = (1 , -a, a 2 , -a 3 ) expt-ay! ) 

02 (Yi) = (1, -p, p 2 , -P 3 ) exp(-pyj) (31) 


Then, the secular determinant that needs to be evaluated at the wall, y = 0, will 
be given by 


f(c) 


0i (0) 0 2 (O) 

D0j (0) D0 2 (O) 


(32) 


B. Results 

As discussed in Section I, there does not exist any mathematical proof on the 
nature of the eigenvalue spectrum for this problem; however, there are two numerical 
investigations on the character of this spectrum. One by Jordinson [6], which is 
considered of a somewhat limited extent, and another more thorough investigation by 
Mack. Since both of these investigations were carried out for the parameters a = 0.308 
and R = 998, these parameters will be used initially for the search procedure to facilitate 
comparison. 

The rectangle in the complex c-plane bounded by 0 < c f < 1.0 and -1.0 <^<0 

was searched for eigenvalues with the parameters a = 0.308 and R = 998. Again, as was 
done for the plane Poiseuille flow, the rectangle was spanned by a number of small 
overlapping circles in such a way as to cover all the desired search area. All of the 
eigenvalues found in this rectangle for this case are shown in Figure 6 and are tabulated 
in Table 3. The value of yj for the initial conditions (31) which was employed in this 
search was 6.0. This represents a length of approximately 1.5 times the boundary layer 
thickness and thus was deemed suitable. 

The structure of the eigenvalue spectrum as shown in Figure 6 is similar to the 
spectra of the plane Poiseuille flow that was discussed earlier. Hence, the classification of 
the spectrum into the A, P, and S families is appropriate for this case also. Note that the 
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Figure 6. All of the eigenvalues found within the indicated boundaries in the c-plane 
for the Blasius boundary layer flow at a = 0.308, R = 998, and yj = 6.0 
(• indicates members of the A family; a indicates members of the 
P family; and A indicates members of the S family). 


spectrum had only one mode possessing c^ > 0, which is the commonly known unstable 

mode. Again, the P and S families occupied regions in the c-plane similar to their 
counterparts for the plane Poiseuille flow with the exception that they did not fall 
strictly on two straight lines. The members of these families seem to lie along two 
neighboring lines that are very close to each other. Hence, the members of the S family 
did not reach a single asymptotic value in c r , rather the two values of c r 0.85 and 

c r « 0.84. However, the modes of this family did extend to as large a value of -c^ as was 
carried out in the search procedure. 
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TABLE 3. EIGENVALUES FOR THE BLASIUS BOUNDARY LAYER AT a = 0.308, 
R = 998, AND y, = 6.0 IN THE REGION OF THE c-PLANE BOUNDED BY 
0.0 < c r < 1.0 AND -1.0 < Cj < 0.0. 


A Family 

P Family 

S Family 

c r 

c i 

c r 

c i 

c r 

c i 

EH 

+0.0079 

0.9913 

-0.0248 

0.8483 

-0.4406 

EH 

-0.1921 

0.9893 

-0.0441 

0.8546 

-0.4852 

EH.. 

-0.2769 

0.9808 

-0.0656 

0.8459 

-0.5357 

0.6863 

-0.3308 

0.9826 

-0.0871 

0.8515 

-0.5764 

0.5572 

-0.3654 

0.9666 

-0.1118 

0.8430 

-0.6372 

0.7939 

-0.4342 

0.9650 

-0.1392 

0.8526 

-0.6754 



0 9501 

-0 1661 

0.8404 

-0.7363 

0.8874 

-0.4148 

vy • y \j a 

0.9474 

-0. 1955 j 



0.9282 

-0.2275 

0.8491 

-0.7823 



0.9229 

-0.2625 

0.8403 

-0.8474 



0.8990 

-0.2978 

0.8483 

-0.8934 



0.8949 

-0.3329 

0.8298 

-0.9611 



0.8726 

-0.3670 





0.8644 

-0.4037 




Thus from the character of the members of the S family it may also be concluded 
that the eigenvalue spectrum of the Orr-Sommerfeld equation for the Blasius boundary 
layer, as posed in expressions (26), (31) and (32) for the numerical integration, 
possesses an infinite and discrete number of eigenvalues. However, it was found that 
the spectrum as shown in Figure 6, and for the parameters used, is not unique. The 
spectrum depended in a strong manner, as far as location and number of eigenvalues 
found, on the value of chosen in the initial conditions (31). However, this 
common character of the spectrum did not affect the overall character, shape, or 
the conclusion on the infinity of eigenmodes of the spectrum. To illustrate this 
dependence of the eigenvalue spectrum on y l5 the spectra at y, = 8.0, 10.0, and 12.0 are 
presented in Figures 7 through 9 with the parameter a and R held fixed. Table 4 lists all 
of the eigenvalues that are shown in Figures 7 through 9 for Cj > -0.7. 
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c. 



Figure 7. All of the eigenvalues found within 0.0 < c r < 1.0 and -0.7 < cj < 0.0 for the 
Blasius boundary layer flow at a = 0.308, R = 998, and y x = 8.0. 


As y r was increased to values greater than 6.0, the number and location of the 
members of the P and S families seemed to vary, while the number and location of the 
members of the A family remained fixed. With the increase of y! , the spacings between 
the modes of the P and S families decreased, resulting in the appearance of new members 
belonging to these families. Furthermore, the juncture point between the P and S families 
shifted to higher values of c r with the increase of y t . This behavior is clearly portrayed 

in Figures 6 through 9 and Tables 3 and 4. These figures also illustrate vividly the shift 
to higher values in c r of the vertical lines along which the modes of the S family are 

located with the increase in y t . It is also clear from these figures that as y t is increased, 
the distance between the adjacent curves, along which the members of the P and S 
families are located, decreased. 
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Figure 8. All of the eigenvalues found within 0.0 < c r < 1.0 and -0.7 < Cj < -0.015 
for the Blasius boundary layer flow at a = 0.308, R = 998, and = 10.0. 


As the modes of the S family shifted to values in c p greater than 0.89, an 

eigenvalue at the location (0.887 - 0.415 i) first appeared. This mode remained in the 
same position for subsequent shifts to the right of the members of the S family. Since 
the position of this mode remained fixed with the increase in y l5 it was concluded that it 
belongs to the A family (in our nomenclature the members of the A family are invariant 
with yj ) and is listed in Table 3 below the dashed line to identify the way it was 
located. Such a result then shows that the lines along which the P and S family lie mask 
any eigenvalue belonging to the A family that lies to the right of these lines; and in order 
to obtain all the members of the A family, the search must be conducted for a very large 
value of yi such that the P and S family lines have shifted to the rightmost location in 
the c-plane. 
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Figure 9. All of the eigenvalues found within 0.0 < c r < 1 .0 and -0.7 < < -0.02 

for the Blasius boundary layer flow at a: = 0.308, R = 998, and y j = 12. 


A comparison of the spectra shown in Figures 6 through 9 with those obtained 
by Jordinson and Mack clearly shows that in both of these investigations only parts of 
the total spectrum were located. By comparing Figure 6 with Figure 1 of Jordinson, it 
appears that he was able to locate some of the modes of the A and P families. However, 
there seems to be no indication that he located any modes belonging to the S family. 
Mack, on the other hand, determined all of the modes of the A family but was not able 
to locate any mode belonging to the P or S families. 
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TABLE 4. EIGENVALUES OF THE P AND S FAMILIES FOR THE BLASIUS 
BOUNDARY LAYER FLOW AT a = 0.308, R = 998, AND y, = 8.0, 10.0, 
AND 12.0 IN THE REGION OF THE c-PLANE BOUNDED BY 
0.0 < c r < 1.0 AND -0.7 < Cj < 0.015. 


yi 

= 8.0 

Yi 

= 10.0 

yi 

- 12.0 

c r 

c i 

c r 

c i 

c r 

c i 

1.00031 

-0.01845 

0.99822 

-0.01645 

0.998503 

-0.01846 

0.99508 

-0.03932 

0.99830 

-0.02194 

0.99873 

-0.02427 

0.98990 

-0.05394 

0.99798 

-0.02979 

0.99754 

-0.03693 

0.99128 

-0.06720 

0.99603 

-0.03780 

0.99529 

-0.04407 

0.98393 

-0.08331 

0.99404 

-0.04833 

0.99556 

-0.05194 

0.98470 

-0.09925 

0.99322 

-0.05679 

0.99441 

-0.06071 | 

0.97726 

-0.11720 

0.99106 

-0.06837 

0.99359 

-0.06935 

0.97623 

-0.13683 

0.99086 

-0.08032 

0.99149 

-0.07907 

0.96817 

-0.15636 

0.98707 

-0.09244 

0.99102 

-0.08889 

0.96599 

-0.17851 

0.98715 

-0.10625 

0.989^3 

-0.09973 

0.95730 

-0.20001 

0.98261 

-0.12046 

0.98818 

-0.11060 

0.95503 

-0.22358 

0.98127 

-0.13572 

0.98562 

-0.12225 

0.94519 

-0.24745 

0.97666 

-0.15078 

0.98500 

-0.13420 

0.94091 

-0.27372 

0.97544 

-0.16706 

0.98208 

-0.14687 

0.93061 

-0.29885 

0.97064 

-0.18387 

0.98109 

-0.15973 

0.92536 

-0.32719 

0.96891 

-0.20154 

0.97836 

-0.17349 

0.91467 

-0.35402 

0.96357 

-0.21970 

0.97652 

-0.18769 

0.90667 

-0.38412 

0.96096 

-0.23902 

0.97334 

-0.20206 

0.88688 a 

-0.40966 3 

0.95572 

-0.25789 

0.97178 

-0.21713 

0.89088 

-0.42565 

0.95271 

-0.27897 

0.96882 

-0.23182 

0.89058 

-0.45890 

0.94642 

-0.29858 

0.96657 

-0.24888 

0.89157 

-0.49188 

0.94341 

-0.32155 

0.96288 

-0.26455 

0.88711 

-0.52771 

0.93603 

-0.34213 

0.96033 

-0.28165 

0.88932 

-0.55986 

0.93373 

-0.36601 

0.95651 

-0.29879 

0.88507 

-0.59833 

0.92495 

-0.39026 

0.95402 

-0.31669 

0.88838 

-0.63201 

0.88736 a 

-0.41476 a 

0.95065 

-0.33398 

0.88364 

-0.67208 

0.92047 

-0.41521 

0.94673 

-0.35375 

0.88714 

-0.70798 

0.91635 

-0.44367 

0.94255 

-0.37232 



0.91409 

-0.47035 

0.93912 

-0.39250 



0.91332 

-0.49612 

0.93521 

-0.41315 
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TABLE 4. (Concluded) 


y i 

= 8.0 

yi 

= 10.0 

yi 

= 12.0 

c r 

■PH 

c r 

c i 

c r 

c i 



0.91084 

-0.52386 

(0.88743 

-0.41472) 



0.91176 

-0.54997 

0.93224 

-0.43476 



0.90944 

-0.57897 

0.93019 

-0.45606 



0.91087 

-0.60619 

0.92848 

-0.47804 



0.90861 

-0.63646 

0.92790 

-0.49909 



0.91011 

-0.66501 

0.92671 

-0.52162 



0.90810 

-0.69627 

0.92669 

-0.54359 



0.90934 

-0.72633 

0.92547 

-0.56691 





0.92584 

-0.58977 





0.92466 

-0.61395 





0.92521 

-0.63772 





0.92406 

-0.66279 





0.92474 

-0.68741 





0.92359 

-0.71338 





0.92436 

-0.73887 


a. These eigenvalues belong to the A family. 


To study the influence of Reynolds number on the character of the eigenvalue 
spectrum, a second case was investigated in which a. was held constant at 0.308 while R 
was increased to 1720. With this variation the qualitative features of the eigenvalue 
spectrum did not change from the previous case; i.e., the spectrum again possessed the A, 
P, and S families. Again, the members of the A family remained invariant with respect to 
y l5 but now the family was comprised of 9 eigenvalues instead of the 7 of the previous 
case. Similarly, the variation of the spectrum with y! followed the same pattern as the 
previous case, as shown in the composite graph of Figure 10. In this figure the eigenvalue 
spectrum for a - 0.308 and R = 1720 is shown for y x = 7.0, 9.0, and 1 1.0. Note the shift 
of the P and S families to the right with the increase in y t . Also note the invariance of 
the location of the members of the A family with respect to yj . Again, the most 
unstable mode belongs to the A family and it is the commonly known unstable mode. 
The eigenvalues that are shown in Figure 10 are listed in Tables 5 and 6. 
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TABLE 5. EIGENVALUES OF THE A FAMILY FOR THE BLASIUS BOUNDARY 
LAYER FLOW AT a = 0.308 AND R = 1720 WITHIN THE REGION OF THE 
c-PLANE BOUNDED BY 0.0 < c r < 1.0 AND -1.0 < cj < 0.0. 


c r 

c i 

0.3360 

0.0041 

0.2412 

-0.2390 

0.4150 

-0.1425 

0.4562 

-0.3187 

0.5775 

-0.2730 

0.6444 

-0.3845 

0.7343 

-0.3571 

0.8189 

-0.4348 

0.8891 

-0.4151 


VI. CONCLUSIONS 


In this report, a purely numerical scheme was presented for locating the temporal 
eigenvalues of the Orr-Sommerfeld equation. The method was especially developed for 
the present investigation to furnish without any ambiguity the information necessary on 
the location of all of the eigenvalues that lie within the investigated region in the 
complex c-plane. The technique was then applied in Sections IV and V to perform the 
desired numerical investigation on the eigenvalue spectrum of the Orr-Sommerfeld 
equation. First the eigenvalue spectra for the plane Poiseuille flow profile were identified 
in Section IV for various values of the Reynolds number while keeping the wave number 
fixed. The search, as expected, did not yield any new information on the nature of the 
eigenvalue spectra for this flow profile. However, the fact that all of the eigenvalues that 
were identified in the search did belong to the spectra and no eigenvalue was found that 
did not belong to it is of greater importance. Hence, it was concluded that the numerics 
of the method had no appreciable influence on the number and location of these 
eigenvalues. Furthermore, extensive tests were conducted to study the effects of the 
numerics on these eigenvalues as they pertain to the integrator step-size, the number of 
orthonormalizations, and the type of the initial value integrator used. The conclusion was 
that all of the eigenvalues that were identified with this technique do belong to these 
spectra. 


Once it was confirmed that the method did yield the correct eigenvalues of the 
problem, the eigenvalue spectrum for the boundary layer profile was investigated in 
Section V. It is concluded from this study that the spectrum obtained from the 
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TABLE 6. EIGENVALUES OF THE P AND S FAMILIES FOR THE BLASIUS 
BOUNDARY LAYER FLOW AT a = 0.308, R = 1720, AND y, = 7.0, 9.0, 
AND 1 1.0 IN THE REGION OF THE C-PLANE BOUNDED BY 
0.0 < c r < 1.0 AND -0.7 < c { < -0.02. 

yf = 7.0 y x = 9.0 yj = 11.00 

Cj. Cj Cj. Cj Cp Cj 

0.99758 -0.01031 0.99919 -0.00766 

0.98937 -0.04736 0.99877 -0.01247 0.99764 -0.02277 

0.99056 -0.05905 0.99699 -0.02390 0.99784 -0.02725 

0.98380 -0.07207 0.99668 -0.02998 0.99662 -0.03287 

0.98380 -0.08568 0.99473 -0.03789 0.99658 -0.03873 

0.97604 -0.10176 0.99442 -0.04448 0.99487 -0.04495 

0.97593 -0.11651 0.99176 -0.05447 0.99502 -0.05128 

0.96734 -0.13360 0.99158 -0.06336 0.99337 -0.05850 

0.96665 -0.15086 0.98887 -0.07241 0.00308 -0.06573 

0.95864 -0.17006 0.98838 -0.08275 0.99053 -0.07366 

0.95524 -0.18942 0.98438 -0.09451 0.99068 -0.08164 

0.94664 -0.20869 0.98397 -0.10576 0.98875 -0.09040 

0.94293 -0.23009 0.97953 -0.11810 0.98814 -0.09910 

0.93305 -0.25114 0.97910 -0.13018 0.98583 -0.10840 

0.92919 -0.27398 0.97471 -0.14325 0.98528 -0.11783 

0.91889 -0.29543 0.97378 -0.15676 0.98223 -0.12785 

0.91371 -0.32048 0.96902 -0.17082 0.98206 -0.13796 

0.90240 -0.34320 0.96744 -0.18540 0.97926 -0.14871 

0.89595 -0.36885 0.96255 -0.20030 0.97831 -0.15981 

0.88319 -0.39172 0.96046 -0.21597 0.97563 -0.17088 

0.97794 a -0.41 420 a 0.95540 -0.23188 0.97439 -0.18257 

0.81941 a -0.43530 a 0.95280 -0.24831 0.97135 -0.19453 

0.87193 -0.43761 0.94761 -0.26491 0.97009 -0.20686 

0.87450 -0.46604 0.94453 -0.28253 0.96681 -0.21930 

0.87013 -0.49752 0.93843 -0.30012 0.96538 -0.23222 

0.87249 -0.52512 0.93518 -0.31846 0.96201 -0.24556 
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TABLE 6. (Concluded) 


y i 

= 7.0 

y i = 

= 9.0 

y i 

= 11.0 

c r 

c i 

c r 

c i 

c r 

c i 

0.86803 

-0.55806 

0.92886 

-0.33681 

0.96011 

-0.25889 

0.87118 

-0.58609 

0.92515 

-0.35657 

0.95682 

-0.27301 

0.86649 

-0.62089 

0.91957 

-0.37475 

0.95439 

-0.28720 

0.87018 

-0.64987 

0.91344 

-0.39687 

0.95073 

-0.30155 

0.86590 

-0.68629 

0.88914 a 

-0.41 507 a 

0.94838 

-0.31631 



0.90667 

-0.41935 

0.94454 

-0.33133 



0.90415 

-0.44353 

0.94187 

-0.34698 



0.90337 

-0.46570 

0.93801 

-0.36246 



0.90100 

-0.48885 

0.93505 

-0.37878 



0.90107 

-0.51036 

0.93100 

-0.39516 



0.89930 

-0.53422 

0.92782 

-0.41256 



0.90069 

-0.55630 

0.88898 a 

-0.41 527 a 



0.89828 

-0.58139 

0.92486 

-0.43022 



0.89951 

-0.60446 





0.89764 

-0.63028 

0.92275 

-0.44833 



0.89908 

-0.65406 

0.92185 

-0.46570 



0.89661 

-0.68109 

0.92012 

-0.48407 





0.92006 

-0.50163 





0.91856 

-0.52051 





0.91901 

-0.53855 





0.91780 

-0.55806 





0.91840 

-0.57669 





0.91724 

-0.59681 





0.91768 

-0.61623 





0.91660 

-0.63685 





0.91730 

-0.65677 





0.91620 

-0.67817 





0.91689 

-0.69866 


a. These eigenvalues belong to the A family. 
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numerical solution of the problem, as was posed, is infinite and discrete. It should be 
remembered, however, that this conclusion is only valid for the numerically approximated 
problem and is not valid for the originally posed problem which is somewhat different. 
This is so because in the process of the numerical approximation to the solution, the 
Orr-Sommerfeld equation was discretised over the finite interval [0,yj], where at the 
point y x the original boundary conditions were replaced with a matched approximated 
analytical solution. Thus due to the restructuring of the problem in the process of the 
numerical solution, the original semi-infinite problem [0,°°) is in effect transformed to a 
finite interval problem. 

This conclusion on the eigenvalue spectrum for the boundary layer profile does 
seem, at first, in contradiction with the conclusion reached by Mack, i.e., that the 
spectrum consisted only of the modes of the A family. But upon a closer examination of 
Mack’s results and taking into account the findings of Jordinson (6) and Benek (9), it 
may be concluded that the eigenvalue spectrum for this problem is as given in Figures 6 
through 10. It should be emphasized here that the latter two results were obtained by 
totally different numerical techniques (Jordinson used the finite-difference technique and 
Benek used the Chebyshev polynomial approximation technique), and thus the modes of 
the P and S families which were found cannot be attributed to numerical inaccuracies 
because in both of these investigations the authors reported finding eigenmodes that 
belong to these families as well as some that belong to the A family. 

In attempting to determine numerically the spectrum for the semi-infinite 
problem, all that can be done is to solve successively a number of finite interval problems 
over the range [0,yj, increasing yj gradually to larger values, and hope to arrive at an 
asymptotic form of the spectrum as y j °°. Such an attempt was undertaken in this 
study and, as far as the numerical calculations allowed, we always found an infinite 
discrete spectrum no matter how large y x was taken. However, as yj became large, the 
number of eigenvalues of the P and S families grew larger and their separation became 
smaller. Also, with the increase of the value of y l5 the lines along which these families 
lie shifted towards the line c r = 1 . This investigation was not able to verify whether or 

not a limit at c r = 1 exists for these families. However, as Figure 1 1 shows, a limit might 
exist for some value of C r <l. 

It may be concluded from these observations that as the interval of numerical 
integration [0,y,] goes to [0,°°), the eigenvalue spectrum for the boundary layer flow 
profile may be composed of two parts: one continuous which is comprised of a straight 
line at c r = 1, and the other discrete and finite comprised of the members of the A 

family. However, for practical purposes in which the problem is solved numerically over a 
finite range of integration, the spectrum will then be considered to be infinite and 
discrete. The only way to resolve this question will be through a rigorous mathematical 
proof. 


34 



1.0 



Figure 1 1. Variation of the value of c r with yi for a typical member of the S 
family for the Blasius boundary layer flow at a = 0.308 and R = 998. 


This investigation may be concluded with some remarks on the utility of the 
numerical search technique used here. A global numerical technique has been outlined for 
use in solving linear differential eigenvalue problems. The technique is general and 
automatic, requiring a minimum of information from the operator. The technique is also 
deterministic and does not require any a priori knowledge on the location of the 
eigenvalues. Furthermore, if high accuracy on the eigenvalues is not desired, the present 
method is sufficient for the solution of the eigenvalue problem. However, when high 
accuracy is required, the present scheme may be used to search for the eigenvalues and, 
once located, their approximate location may be input into a local method that will 
converge fast enough to the desired accuracy on the eigenvalues and the eigenfunctions. 
Thus, this method could save considerable time by eliminating the trial and error 
procedure and the guesswork that might be involved when using a local technique. 
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